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The hydrodynamic core of the Terminal Area Simulation System (TASS) is evaluated 
against different benchmark cases. In the absence of closed form solutions for the equations 
governing atmospheric flows, the models are usually evaluated against idealized test cases. 
Over the years, various authors have suggested a suite of these idealized cases which have 
become standards for testing and evaluating the dynamics and thermodynamics of 
atmospheric flow models. In this paper, simulations of three such cases are described. In 
addition, the TASS model is evaluated against a test case that uses an exact solution of the 
Navier-Stokes equations. The TASS results are compared against previously reported 
simulations of these benchmark cases in the literature. It is demonstrated that the TASS 
model is highly accurate, stable and robust. 


I. Introduction 

T he Terminal Area Simulation System (TASS) is a state-of-the-art, cloud-resolving, large eddy simulation 
atmospheric flow model. TASS computes the primitive variables non-hydrostatic equations in three dimensions 
and the model is capable of resolving flows at multiple spatial (spatial resolutions varying from less than a meter to 
2km in the horizontal) and temporal scales (few minutes in case of wake decay simulations to hours for studying 
convective phenomena). The model also solves prognostic equations for potential temperature, water vapor, cloud 
droplets, ice crystals, rain, snow and hail. Phase changes between liquid water and ice as well as submodels for 
cloud and precipitation microphysics are included. Subgrid scale diffusion is parameterized via Smagorinsky-type 
turbulence closure (Smagorinsky 1963) and surface layer processes are computed based on Monin-Obukhov 
similarity theory (Stull 1997). 

The TASS model has been extensively validated in the past against experimental data and applied to a diverse 
set of problems ranging from the simulation of shallow cumulus to supercell storms, including convective 
phenomena such as downbursts, gust fronts and hail storms (Proctor 1987; Proctor and Bowles 1992; Proctor et al. 
2002). Although, TASS was initially developed to simulate convective processes, it has been used successfully to 
simulate the transport and dispersion of hazardous materials from nuclear weapon bursts (Bacon and Sarma 1991). 
TASS simulations of large urban fires have been used to define source terms for mesoscale models (Bacon et al. 
1986). More recently, TASS has been employed to study the decay and transport of wake vortices generated by 
various aircraft under different conditions of atmospheric stability and turbulence (Proctor 1996; Han et al. 2000a; 
Han et al. 2000b; Proctor 2009). 

In this study the dynamical core of the TASS model is evaluated against benchmark cases and compared with 
results from other numerical models such as the National Center for Atmospheric Research’s (NCAR) state-of-the- 
art Weather Research and Forecast (WRF) model and the Advanced Regional Prediction System (ARPS). A 
detailed description of the WRF model is given in Klemp et al. (2000) and the ARPS model is described in Xue et 
al. (2000). Over the years a suite of idealized test cases have been developed for the validation and verification of 
atmospheric models. Although these cases can evaluate a model mostly in qualitative terms, they nonetheless 
provide a useful measure to assess a model’s capability to simulate the fundamentals of atmospheric hydrodynamics 
and thermodynamics. Three idealized test cases were used in this study - the non-linear density current simulation 
(Straka et al. 1993); convection of a warm bubble in a neutral atmosphere (Bryan and Fritsch 2002); and the 
propagation of non-hydro static inertia gravity waves (Skamarock and Klemp 1994). In addition to the idealized test 
cases mentioned above, the Beltrami flow problem was also computed with TASS and compared against the exact 
solution (Shapiro 1993). 
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II. Governing Equations 

The complete governing equations of the TASS model consist of conservation of mass, momentum and energy 
in three dimensions that include the subgrid-scale diffusion processes and comprehensive cloud and ice 
microphysics (Proctor 1987). In this study, the dynamical core of the TASS model is evaluated and therefore, the 
diffusion processes are ignored and the atmosphere is assumed to be dry. The simplified equations can then be 
written as follows: 


a., + 7 / y a, +g(H _ l)s 

dt p 0 dx i dxj dxj 

(1) 

where, u x and wj are the velocity components, x x and Xj are the spatial dimensions, p is the pressure, g is the 
acceleration due to gravity, and 8 is the Kronecker delta. H is the density ratio defined by: 

/ \C V /C p 

H _ A> _ 0 p 0 
P #0 l P J 

(2) 

In Eq. (2), C p (= 1004 J K" 1 kg" 1 ) and C v (= 111 J K' 1 kg' 1 ) are the specific heats of air at constant pressure and 
volume respectively. p 0 is the base-state pressure, p 0 is the base-state density and 6 {) is the base-state potential 
temperature. The horizontal velocities, atmospheric pressure, potential temperature and density are decomposed 
into initial time- invariant hydrostatic base-states (u 0 ,v o ,p o ,0 o ,p o ) which are a function of height only and the time- 
variant perturbations (u\v\ p\6\ p r ) which are a function of both space and time: 

u(x, y,z,t) = u 0 (z) + u '(x, y, z, t) 

(3) 

v(x, y,z,t) = v 0 (z) + v'(x, y,z,t ) 

(4) 

p(x, y, z, t) = p 0 (z) + p'(x, y, z, t) 

(5) 

6{x, y, z, t) = 0 o (z) + d'(x, y, z, t) 

(6) 

p(x, y, z, t) = p 0 (z) + p\x, y, z, t ) 

(7) 

The prognostic equation for pressure deviation is approximated by: 


dp' C duj 

ir+c: p *;- p ‘ su ^ 

(8) 


The conservation of energy is given in terms of potential temperature with the assumption of a dry adiabatic 
atmosphere: 


dO du .6 du. 

— = J — + 6 ^~ 

dt dxj dxj 

The system is closed by an equation of state for pressure: 

p = r ,jp t 

In Eq. (10), T is the air temperature and R d (= 287 J K' 1 kg' 1 ) is the gas constant for dry air. 


(9) 


( 10 ) 
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III. Numerical Methods 

The numerical methods used in TASS are described in detail by Proctor (1987; Proctor 1996). A brief 
description is given here for the sake of completeness. The equations are discretized using a quadratic-conservative, 
fourth-order finite-differences in space for the calculation of momentum and pressure fields (Proctor 1996) and the 
third-order Leonard scheme (Leonard 1995) is used to calculate the transport of potential temperature. A Monotone 
Upstream-centered Scheme for Conservation Laws (MUSCL)-type scheme after van Leer (van Leer 1979; Ahmad 
and Proctor 2011) has also been implemented in TASS for the transport of scalar quantities. The TASS 
computational mesh uses the Arakawa-C grid staggering for specifying velocities and thermodynamic quantities 
(Arakawa and Lamb 1977). The Klemp-Wilhelmson time-splitting scheme (Klemp and Wilhelmson 1978) is used 
for computational efficiency in which, the higher- frequency terms given in the left hand side of Eq. (1) and (8) are 
integrated by enforcing the CFL criteria to take into account sound wave propagation due to compressibility effects. 
The remaining terms in Eq. (2) and (8) and Eq. (9) are integrated using a larger time step that is appropriate for 
anelastic and incompressible flows. Non-reflecting Orlanski boundary conditions (Orlanski 1976) are imposed at 
the outflow boundaries. A sixth-order filter is used to damp-out spurious oscillations that may arise due to the use of 
centered-differencing of momentum and pressure terms. The model code is written in FORTRAN 77 and has been 
fully parallelized for distributed computing platforms using the Message Passing Interface (MPI). 

IV. Results 

In this section, the different benchmarks cases simulated to evaluate the hydrodynamic core of TASS are 
described in detail. 

A. Straka Density Current 

The non-linear density current benchmark (Straka et al. 1993) is presented in this section. The domain in this 
case was bounded from -20km to 20km in the x-direction; from 0 to 5km in the y-direction; and from 0 to 6.4km in 
the z-direction. The mesh resolution was set to 50m in both the x- and z-directions and 1km in y-direction. Outflow 
boundary conditions were used in the x-direction and periodic boundary conditions were used in the y-direction. 
The domain was initialized for a neutral atmosphere by setting the potential temperature at 300 K. The initial u- 
velocity and the w- velocity were set to zero. The temperature profile was defined using the following relation: 


r(z) = 300-^ 

C 


(ii) 


A cold bubble was initialized by adding a perturbation in the temperature field using the following relation: 



fo.o 


AT = < 

1 

1-15.0 

cos(^L)+1.0 


l 

2 


if L> 1.0, 
if L <1.0 


( 12 ) 


where, L = ^ j(x - x c / x r ) 2 + (z - z c / z r ) 2 , x c = 0, z c = 3km, x r = 4km and z r = 2km. As prescribed in Straka (1993), a 

constant eddy viscosity/diffusivity ( K m = K h = K = 75 m 2 /s) was assumed for both the momentum and potential 
temperature fields. The simulation was run for time = 900s. The potential temperature perturbation (K), ^-velocity 
(m/s) and the w- velocity (m/s) fields at time = 900s are shown in Figure 1. The figures show a portion of the actual 
computational domain from approximately 0 to 20km in the x-direction and 0 to 6000m in the z-direction. The 
simulation results shown in Figure 1 are in good agreement with the Straka benchmark, both qualitatively and 
quantitatively (see Table 1). In Table 1, REFC refers to the fully-compressible reference solution, REFS is the fully 
compressible reference solution on a staggered grid and REFQ is the reference solution using a quasi-compressible 
model. All reference solutions (REFC, REFS and REFQ) reported in Straka et al. (1993), were on grids with 
uniform resolution of 25m. /-waves refers to a Godunov-type scheme reported in Ahmad and Lindeman (2007). 
The data in Table 1 are for the right half of the computational domain (shown in Figure 1). The location of the front 
for the TASS solution is in terms of potential temperature and is given by the x-coordinate of the cell in the lowest 
level with a potential temperature of 299.99 K. TASS solution also compares well with other reported simulations 
of the benchmark (e.g., see Janjic et al. 2001; Xue et al. 2000; Giraldo and Restelli 2008). 
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Table 1: Comparison of the TASS results with the Straka Density Current Benchmarks 


Variable 

TASS 

Ax = Az = 50m 

f-waves 
Ax = A z — 50m 

REFC 

Ax = A z — 25m 

REFS 

Ax = Az = 25m 

REFQ 

Ax = A z — 25m 

p' max (mb) 

1.31 

1.26 

2.87 

2.49 

1.74 

p' min (mb) 

-6.02 

-6.27 

-5.14 

-5.55 

-5.21 

0 'max (K) 

0.0359 

8.92E-03 

0.00 

0.00 

0.00 

0' ml „(K) 

-9.24 

-9.82 

-9.77 

-9.77 

-10.00 

«max (m/s) 

37.35 

34.44 

36.46 

35.02 

34.72 

^min (m/s) 

-14.49 

-15.74 

-15.19 

-16.32 

-15.31 

H’max (m/s) 

12.17 

13.62 

12.93 

13.28 

13.04 

W m in (m/s) 

-15.16 

-16.36 

-15.95 

-16.11 

-16.89 

front location (m) 

15596.3 

15525.00 

15537.44 

N/A 

15509.17 
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Figure 1. Non-Linear Density Current Benchmark. Potential temperature perturbation (K) - top, w-velocity 
(m/s) - middle, and w-velocity (m/s) - bottom. Time = 900s. Av=50m and Az=50m. The maximum and 
minimum values of each field are given on the top-left corner of each figure. Dashed lines in the velocity 
graphics indicate negative values. Only the right half of the domain is shown in the figure. 
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B. Convection in Neutral Atmosphere 

There are many variants of this test case, e.g., see Koracin et al. (1998), Ahmad and Lindeman (2007) and Bryan 
and Fritsch (2002). TASS was evaluated using the setup similar to that described by Bryan and Fritsch (2002). The 
domain for this case was bounded from 0 to 20km in the x-direction; from 0 to 5km in the y-direction and from 0 to 
10km in the z-direction. The mesh resolution was set to 100m in the x- and z-directions and 1km in the y-direction. 
Outflow boundary conditions were used in the x-direction and periodic boundary conditions were used in the y- 
direction. The top and bottom boundaries were set to solid walls (Bryan and Fritsch set wall conditions on the 
lateral boundaries as well). The domain was initialized for a neutral atmosphere at = 300 K in hydrostatic 
balance. The initial w-velocity, v- velocity and the w - velocity were set to zero. A potential temperature perturbation 
in the form of a smooth cosine function was initialized near the lower boundary of the domain: 


0\x,z,t = 0) = AO 0 cos 2 



where, 


L = 



+ 


/ \2 
Z ~ Z c 

K Z r , 


(13) 


(14) 


The potential temperature perturbation maximum A# 0 , was set to 2K. The center of the thermal perturbation was 
prescribed at x c = 10km and z c = 2km. The radius of the perturbation was set to 2km (i.e., x r = z r = 2km). The model 
was run for time = 1000s. The introduction of a warm potential temperature perturbation generates strong updrafts 
in the thermal core accompanied by downdrafts on either sides of the thermal. This is a highly non-linear case and 
as the thermal rises, large gradients in potential temperature develop at the top of the warm bubble. The TASS 
simulation results of potential temperature perturbation after time = 1000s is shown in Figure 2 along with results 
from the simulation by Bryan and Fritsch (2002). The TASS computed vertical velocity field after time = 1000s is 
shown in Figure 3 along with results from the simulation by Bryan and Fritsch (2002). TASS was able to accurately 
simulate the convection in neutral atmosphere. The symmetry is well maintained in both the potential temperature 
and velocity fields and the results compare well qualitatively as well as quantitatively with the data reported by 
Bryan and Fritsch (2002). Bryan and Fritsch used a fifth-order upwind-biased scheme which may explain the slight 
differences in their solution. 



x (km) 



Figure 2. Convection in Neutral Atmosphere. Potential Temperature (K) field is shown in the figure. 
Comparison of the TASS computed solution (left) with the fifth-order upwind scheme (right) reported in 
Bryan and Fritsch (2002). Time = 1000s. Ax=100m and Az=100m. 
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Figure 3. Convection in Neutral Atmosphere. Vertical velocity (m/s) field is shown in the figure. 
Comparison of the TASS computed solution (left) with the fifth-order upwind-biased scheme (right) reported 
in Bryan and Fritsch (2002). Time = 1000s. Ax=100m and Az=100m. 


C. Non-hydrostatic Inertia Gravity Waves 

The Skamarock-Klemp (1994) test for simulating inertia-gravity waves on the non-hydrostatic scale is described 
in this section. The domain was bounded from 0 to 300km in the x-direction, from 0 to 5km in the y-direction and 
from 0 to 10km in the z-direction. Periodic boundary conditions were used in the lateral and the top and bottom 
boundaries were set to solid walls. The potential temperature in the vertical was initialized by a constant Brunt- 
Vaisala frequency N= 10" 2 s -1 . The waves were excited by an initial potential temperature perturbation given by: 


#(x,z,0) = A0 O 


sin(;rz/ H ) 
l + (x- x c ) 2 / a 2 


(15) 


The amplitude of the initial potential temperature perturbation, A0 o was set to 10" 2 K. The height H of the 
domain was 10km, the perturbation half width was a = 5km. The perturbation was initialized at x c = x max /3, where, 
x max is the width of the domain (300km). The initial w-velocity and w-velocity were set to zero. The model was run 
for time = 3000s. The mesh resolution in this simulation was set to 1km in the horizontal and 100m in the vertical. 
Due to the presence of gravity waves, the large time step used for advecting potential temperature in the time- 
splitting scheme was set to 2.5 times the acoustic time step. 

Figure 4 shows the comparison of the centerline data at z = z max /2 for x = [0:300km] for TASS model, the WRF 
fifth-order upwind scheme and the WRF second-order centered scheme (see Ahmad and Lindeman 2007; Ahmad et 
al. 2007 for details on the WRF simulations). The TASS profile matches the WRF fifth-order solution quite well, 
whereas large dispersion and phase errors can be seen in the WRF second-order solution. The TASS solution on the 
other hand exhibits no noticeable diffusion, phase or dispersion errors. The potential temperature perturbation field 
for the TASS model after 3000s into the simulation is shown in Figure 5. 

It should be noted that the solution in Skamarock and Klemp (1994) assumes a Boussinesq atmosphere and the 
computed results presented here are obtained from compressible flow models. The TASS results also compare well 
with other reported solutions of this benchmark (Ahmad et al. 2007; Ahmad and Lindeman 2007; Girlado and 
Restelli 2008). 
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Figure 4. Non-hydrostatic Inertia-Gravity Waves. Comparison of the TASS model with the WRF fifth-order 
upwind scheme and the second-order centered scheme. Theta perturbation (K) values between x = 0 & 
300km, for z = 5km are shown in the figure. 300x100 cells (Av=lkm and Az=100m). Time = 3000s. 
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Figure 5. Non-hydrostatic Inertia-Gravity Waves. TASS computed potential temperature perturbation (K) 
field. Time = 3000s. 300x100 cells (Ac=lkm and Az=100m). 
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D. Viscously Decaying Beltrami Flow 

A few closed-form solutions to the atmospheric flow equation sets have been found to exist under special 
circumstances. One such test case of the Beltrami flow was suggested by Shapiro (1993) for which an exact 
analytical solution exists. This test case is invaluable in verifying quantitatively, the accuracy and stability of the 
spatial and temporal discretization schemes implemented in an atmospheric flow model. The exact solution of the 
Beltrami flow is described in detail by Lilly (1986), Shapiro (1993), and Wang (1990). The simulations of the 
Beltrami flow using an earlier version of TASS are given in Switzer (1996). Please note that all of the previous 
cases presented in this paper described phenomena in two dimensions whereas the Beltrami test case is in three 
dimensions. The computational domain in this case comprised of 89 x 59 x 44 cells and the mesh resolution was set 
to 3m in the x-direction, 2m in the y-direction and lm in the z-direction. The time-varying viscous Beltrami flow 
field is given by: 


u(x,y,z,t) ■ 


k 2 +i 


- \Al cos(Ax) sin(/y) sin(mz) + mk sin(Ax) cos(/y) cos(mz)]exp(-vT 2 ^) . 


(16) 


v(x, y, z, t ) = ^ \A k sin (kx) cos (ly) sin (mz) - ml cos (kx) sin (ly) cos(raz)]exp(-v2 2 ^) , 


k 2 +/ 2 


(17) 


w(x,y,z,t ) = Acos(kx)cos(ly)sin(mz)exp(-vA 2 t ) , 


(18) 


p(x,y,z,t) = p s 



u 2 +v 2 + w 2 
2 


+ gz , 


where, 


A 2 =k 2 + 1 2 +m 2 . 


and, 


k 





(19) 


( 20 ) 


( 21 ) 


In Eq. (16)-(21), k , / and m are the wave numbers. p s is the stagnation pressure at ground level and A is the 
amplitude of the vertical velocity. The gradients of u-, v- and w-velocity can be increased or decreased by changing 
the wave numbers k , /, and m. vis the viscosity, g is the acceleration due to gravity and p is the fluid density. 
Initial conditions for the three components of the velocity field (u, v, and w) and the pressure, p can be obtained by 
setting the time, t = 0. For the TASS simulations presented in this section, A = 2m/s, Z x = 267m, L y = 1 18m and L z = 
44m are the domain bounds in the x-, y- and the z-direction respectively. Although, the exact Beltrami solution 
assumes an incompressible fluid, its use for validation of compressible flow models can be justified if the height of 
the computational domain is set such that it is much smaller than the density scale height (Shapiro 1993). The use of 
periodic boundary conditions in all directions is required for this test case. Two simulations were run with different 
values of viscosities (v=lm 2 /s and v=0.1m 2 /s). Periodic boundary conditions were used in all directions and the 
final integration time was set to 41s as in Shapiro (1993). 

Figure 6 shows the RMS error in domain maximum vertical velocity and Figure 7 shows the RMS error in 
domain total kinetic energy as a function of time. The errors for both the decaying and the weakly decaying cases 
are minimal. Figure 8 shows the computed u-, v-, and w-velocity fields in the x-y plane at z = 10.5m. The difference 
between the exact and the computed solutions is shown in Figure 9. It can be seen that the errors in the TASS 
solution are very small and the fidelity of the solution is comparatively better than the Advanced Regional 
Prediction modeling System (ARPS) solution given in Shapiro (1993). Figure 10 shows the i/-velocity and w- 
velocity values between x = 0 and 267m, fory = 60m and z = 10.5m. The computed TASS results are in excellent 
agreement with the exact solution. 
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Figure 6. Viscously Decaying Beltrami Flow. Time history of RMS error in domain maximum w-velocity 
(top) and domain total kinetic energy (bottom) using different values of viscosities - the decaying case (v = 
lm 2 /s) and the weakly-decaying case (v= 0.1m 2 /s). 89x59x44 cells. Av=3m, Ay = 2m, and Az=l m. 



Figure 7. Viscously Decaying Beltrami Flow. Time history of RMS error in domain maximum w-velocity 
(top) and domain total kinetic energy (bottom) using different values of viscosities - the decaying case (v = 
lm 2 /s) and the weakly-decaying case (v= 0.1m 2 /s). 89x59x44 cells. Ar=3m, Ay = 2m, and Az=lm. 
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Figure 8. Viscously Decaying Beltrami Flow. TASS solution: w-velocity in m/s (top panel); v-velocity in m/s 
(middle panel) and the >c-velocity in m/s (bottom panel). The x-y plane is shown at z = 10.5m. Viscosity 
coefficient, v = lm 2 /s. Time = 41.02s. 89x59x44 cells. Av=3m, Ay = 2m, and Az=lm. 
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Figure 9. Viscously Decaying Beltrami Flow. Difference between computed and exact solutions: w-velocity 
in m/s (top panel); v-velocity in m/s (middle panel) and the w-velocity in m/s (bottom panel). The x-y plane is 
shown at z = 10.5m. Viscosity coefficient, v = lm 2 /s. Time = 41.02s. 89x59x44 cells. Ar=3m, Ay = 2m, and 
Az=lm. 
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Figure 10. Viscously Decaying Beltrami Flow. Comparison of computed and exact solutions, w-velocity (m/s) 
values between x = 0 & 267m, for y = 60m and z = 10.5m are shown in the top panel, w-velocity (m/s) values 
between x = 0 & 267m, for y = 60m and z = 10.5m are shown in the bottom panel. Viscosity coefficient, v = 
lm 2 /s. Time = 41.02s. 89x59x44 cells. A*=3m, Aj = 2m, and Az=l m. 


12 

American Institute of Aeronautics and Astronautics 


V. Summary 

The TASS hydrodynamics was evaluated against three benchmark cases in two dimensions and against an 
analytical test case in three dimensions with excellent results. The higher-order spatial and temporal discretization 
schemes implemented in TASS exhibit minimal dispersion, phase and diffusion errors. The simulation results also 
compare well with other models. In the inertia-gravity wave case, for example, the TASS model performs at par 
with the fifth-order upwind biased scheme implemented in WRF and out-performs WRF’s second-order centered 
scheme which exhibits large dispersion and phase errors. In the case of Beltrami flow, again the TASS model 
performed exceptionally. In terms of accuracy, TASS also out-performs the ARPS model for Beltrami simulation 
(see Shapiro 1993). Due to its higher order spatial accuracy and minimal dispersion and phase errors, TASS is well 
suited for simulating complex turbulent flows and convection dominated flows. 
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